Eur ophysics Letters PREPRINT 



The rich behavior of the Boltzmann equation for dissipa- 
tive gases 

M.H. Ernst 1 , E. Trizac 2 and A. Barrat 3 

1 Theoretische Fysica, Universiteit Utrecht, Postbus 80.195, 3508 TD Utrecht (The 
Netherlands) 

2 Theoretical Biological Physics, UC San Diego, La Jolla CA 92093 (USA ) and LPTMS 
(UMR CNRS 8626), Univer site Paris- Sud, 91405 Orsay (France) 

3 LPT (UMR CNRS 8627), Universite Paris-Sud, 91405 Orsay (France) 



PACS. 45.70.-n - Granular systems. 
PACS. 05.20.Dd - Kinetic theory. 

PACS. 81.05.Rm - Porous material; granular materials. 



Abstract. - Within the framework of the homogeneous non-linear Boltzmann equation, 
we present a new analytic method, without the intrinsic limitations of existing methods, for 
obtaining asymptotic solutions. This method permits extension of existing results for Maxwell 
molecules and hard spheres to large classes of particle interactions, from very hard spheres to 
softer than Maxwell molecules, as well as to more general forcing mechanisms, beyond free 
cooling and white noise driving. By combining this method with numerical solutions, obtained 
from the Direct Simulation Monte Carlo (DSMC) method, we study a broad class of models 
relevant for the dynamics of dissipative fluids, including granular gases. We establish a criterion 
connecting the stability of the non-equilibrium steady state to an exponentially bound form 
for the velocity distribution F, which varies depending on the forcing mechanism. Power laws 
arise in marginal stability cases, of which several new cases are reported. Our results provide 
a minimal framework for interpreting large classes of experiments on driven granular gases. 



Granular materials are ubiquitous in natural phenomena and widely used in industrial pro- 
cesses. Yet, a fundamental understanding of their properties still challenges scientists [1,2]. 
Although a large number of particles is involved, a thermodynamic-like description remains 
elusive and has been thwarted by energy losses through internal dissipation, and energy sup- 
plies from non-thermal sources. Granular media are therefore intrinsically far from equilib- 
rium, and a paradigm for open dissipative systems. One can distinguish two types of granular 
fluids [1]: quasi-static flows (solid-like), where gravity, friction and stress distributions control 
the dynamics, and rapid granular flows (gas- like), where the dynamics is mostly governed by 
ballistic motion in between inelastic collisions. This letter deals with granular gases. 

Even dilute granular gases differ from their molecular counterparts to a significant extent [1, 
2], and their particle velocity distribution is emblematic of the difficulties of a statistical 
description. The experimentally measured velocity distribution F(v) generically deviates 
from the Maxwell-Boltzmann behavior, and its functional form -often fitted by stretched 
exponentials- depends on material property, on the specific geometry considered and on the 
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forcing mechanism which compensates the collisional loss of energy [3-7]. A similar picture 
emerges from numerical simulations and analytical studies [5,8-25], where in addition to 
the more common stretched exponential behavior, power-law distributions have also been 
reported [16,17,25]. This rather vast but disparate body of knowledge calls for rationalization. 

Our goal is to propose a general and unified framework for interpreting experimental 
results, as well as to generalize existing theoretical results in two new and important directions: 
(i) to more general repulsive interactions, parametrized by an exponent v , which include the 
known results for Maxwell molecules [y = 0) and hard spheres [y = 1), and (ii) to more 
general forcing mechanisms, such as nonlinear negative friction forces, parametrized by an 
exponent 6, including free cooling (9 = 1), as well as the white noise driving force. In 
addition, since the Maxwell-Boltzmann distribution seems to have no analog in steady states 
of dissipative gases, we want to identify the generic trends for the velocity distribution and 
study its properties in this framework. To this aim, we develop a global and quantitative 
cartoon of the different effects at work, that lead, upon changing certain key parameters, to 
the wide range of behaviors alluded to above. 

In this letter, we develop a new analytic method for determining deviations from Gaus- 
sian behavior, that is free of the intrinsic limitations of previous approaches. For instance, 
the known Fourier transform method for obtaining high energy power law tails [16, 17, 27] 
is intrinsically restricted to Maxwell molecules, where the nonlinear collision term is a con- 
volution product. Moreover, the method for obtaining stretched exponential tails in [11,17] 
was based on a one-stroke asymptotic approximation, whereas our new method provides a 
systematic approximation scheme yielding large sub-leading corrections. We also show that 
the large velocity tail of F(v) is characterized by an exponent b, which governs the stability 
of the non-equilibrium steady state. When this state is a stable fixed point of the dynamics, 
b > and F is of stretched exponential form ~ exp[— v b \. We also emphasize the relevance 
of subdominant corrections that must be taken into account for a comparison of numerical 
or experimental data with the theory. For cases of marginal stability where b vanishes, new 
classes of power law tails v~ s appear. An immediate experimental consequence is that such 
power law tails will be elusive, since they can only arise when material properties are fine 
tuned. Small variations will either drive the system into an unstable state, or into a stable 
state with a stretched exponential tail. 

We use the simplest model for rapid granular flows, i.e. the nonlinear Boltzmann equa- 
tion, that corresponds to a regularly 'randomized' gas (by shaking, fluidization, vibrational 
techniques, or bulk driving through alternating electric fields). The dynamics is described as 
a succession of uncorrelated inelastic binary collisions, modelled by soft spheres with collision 
frequency, g<,(g,d) ~ g u \% ■ n\ a , and a coefficient of normal restitution a < 1. The collision 
law (vi, V2) — > (v[, v 2 ) reads: = vi — p(g-n)n where g = vi — V2, p = 1 — q = \(l + a) and 
n is a unit vector parallel to the impact direction connecting particles 1 and 2. The expres- 
sion for v' 2 follows from momentum conservation. Momentum transfer at a distance can also 
take place through long range interactions [4,5]. Collisions dissipate kinetic energy at a rate 
cc i(l —a 2 ) = 2pq. The quantity q(g 1 -d) is the differential scattering cross section, where expo- 
nent a models the angular dependence, and g = g/g is a unit vector, parallel to g. For elastic 
particles interacting via a soft sphere potential U(r) oc r~ a , one has v = 1 — 2(d— l)/a, where 
d is the space dimension. The exponents v = a = 1 model ordinary hard-sphere behavior 
(a — > 00), and v — defines so-called Maxwell models. The hallmark of Maxwell interactions 
is a collision frequency independent of the relative velocity of the colliding pair. In Ref. [5] 
this property has been found experimentally, and explained theoretically for experiments on 
non-magnetic particles in a viscous fluid, and magnetic grains with dipolar interactions in air 
in a 2-D geometry where U(r) = r~ 2 . Because of this hallmark property, Maxwell models 



M.H. Ernst, E. Trizac and A. Barrat: The rich behavior of the Boltzmann equation for dissipative gases3 



are particularly convenient for analytical purposes. Note also that the importance of material 
properties has barely been taken into account in previous studies. Here we do so in a nec- 
essarily schematic but physically relevant manner through the parameters a and a, that 
encode material properties. 

The study of the homogeneous systems is not only a useful starting point, it is also relevant 
to experiments with bulk driving [5,6], and a significant fraction of the experiments [3-7] has 
been carried out on spatially uniform systems. The homogeneous non-linear Boltzmann reads 

d t F(v,t)+FF = J dv 1 dv 2 dn g ^|g-n|' T F(v 1 ,t)F(v 2 ,t)[ ( 5(v-v' 1 )- ( 5(v-v 1 )], (1) 

where the r.h.s defines the collision kernel I(v\F) with the usual gain/loss structure and Dirac 
functions to enforce the collision law; J indicates that the n -integration is an average over 
the pre-collision hemisphere, g n < 0, with a uniform weight normalized to unity. The forcing 
term TF represents the energy supply. We consider here the generic form 

TF = 7 <9 V • (vv e - 1 F) - DdlF (2) 

with either (7 = 0, D > 0) or (7 > 0, D — 0), corresponding respectively to stochastic White 
Noise (WN), or deterministic nonlinear Negative Friction (NF), where 8 > is a continuous 
exponent. The change in velocity due to the latter force is proportional to tr — 1 v. For instance, 
for 9 > 1 it injects selectively energy in the high energy tail of F(v). A physical realization of 
NF with 9 > 1 is provided by an electrically forced gas with a charge distribution. Energy is 
predominantly injected in the particles with a high charge, which in turn reach larger typical 
velocities. Recently a somewhat different device with a similar property has been proposed 
in Ref. [25]. The special case 9 — 1 corresponds to the Gaussian thermostat [12], allowing 
to study the long time scaling regime of an unforced system (so-called free cooling), while 
9 = models gliding friction, also called gravity thermostat [12]. It corresponds to a force of 
constant strength acting in the direction of the velocity. WN forcing has been frequently used 
in analytical and numerical studies [11, 12, 14, 19,23], and more importantly, was shown to be 
experimentally relevant for two dimensional experiments with energy injected from a vibrating 
rough plate [6,7] . In such cases individual grains undergo random changes in direction, that are 
mimicked by white noise forcing. Existing experimental evidence also shows that stochastic 
forcing is frequently a valuable first approximation [4, 14]. We therefore have a restricted 
'phase space' to explore: (v, a, a) for WN and (v, a, a, 9) for NF. Here 7 and D are irrelevant 
constants that set the energy scale. 

While the dynamics defined by Q and © always admits a non-equilibrium steady state, 
the first question to address is the stability of this steady state. Our analysis assumes a rapid 
approach of F to a scaling form F(v,t) = (l/vo(t)) d f(v/vo(t)) (for proofs, see [18]), where 
the r.m.s. velocity v is defined through the granular temperature, J dvv 2 F(v,t) = ^dvf^t), 
and Eq.QJ implies: 

dvl/dt oc 1 - (v (t) / v (oo)) 2b (b WN = 1 + v/2) 
dvo/dt ocug(t) [l - (v (t)/v {oo)) b ] (b NF = 1 + v - 9) (3) 



The solution v therefore approaches the fixed point Wo(oo) if b > 0, and conversely moves 
away from vq(oo) if b < 0, in which case the steady state is an unstable fixed point of the 
dynamics. 
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Fig. 1 - Comparison of the velocity distribution obtained from Monte Carlo simulations with the 
asymptotic predictions. Left: WN forcing for d = 2, a — 1/2, a = v = —1/2, where b = 3/4, b' = 
9/16, x = 0. DSMC data are plotted as /(c) (dashed-dotted line) and exp [-P'c b ]/(c) (solid line) vs 
x — f3c b , and compared with the theoretical prediction e~ x (dashed line). Right: NF forcing with 
d = 1, a = 0, v = 1/2. The DSMC data for /(c) vs c b are compared at various 8 with the prediction 
exp[— fic ] (dashed lines). f3,P' have been measured in the DSMC simulations from their definition 
0. Fig. 2: WN-driven system at stability threshold (d = 2, v = —2, a = 1) for different values of the 
restitution coefficient a. The dashed lines correspond to the theoretically predicted power-law tails 
l/c a . The inset shows the algebraic exponent obtained through the numerical solutions of Ea. Ill 111 . 

In stable steady states the integral equations for /(c) can be derived from JIJ, an d read 
for WN and NF driving respectively: 

mf) = -J^sdlf = -^ d ((9 v+2 ))dlf 

W) = (M^ c .(cc*/) = ^$^d c -(cc°f) (4) 

The phenomenological constants D and 7 have been eliminated with the help of the steady 
state conditions in The averages ((■ ■ •)) and (■ • •) are respectively taken over the unknown 
/(ci)/(c 2 ) and /(c), and g = \ci — c 2 |. Moreover A 2 = 2pq/3 a+2 and (3 S = f dn |a • n\ s / J dn 
is an average of | cost9| s over a <i-dimensional solid angle. 

Next we discuss the new criterion relating stretching exponents to stability. From the one 
stroke asymptotic approximation of [11] an interesting feature already emerges: in the stability 
region (b > 0), /(c) is of stretched exponential form with an a priori known exponent b from 
(PJ, while in cases of marginal stability (b — > + ), / ~ c~ a is of power law type with an a 
priori unknown exponent a. As expected, b decreases with v, since a tail particle with velocity 
c> 1 suffers collisions at a rate c v ; the slower the rate, the slower the particle redistributes 
its energy over the thermal range c ~ 1, which results in an increasingly overpopulated high 
energy tail. When v is further decreased so that b changes sign, the tail is no longer able 
to equilibrate with the thermal "bulk" , and the system cannot sustain a steady state. A 
similar intuitive picture may be developed with respect to 9 in the NF cases. While 9 = 1 is 
equivalent to a linear rescaling of velocities, a nonlinear rescaling with 9 > 1 selectively puts 
more energy in the tail particles, thus leading to an increased overpopulation, as b decreases 
(see ©)■ The reverse scenario applies for 9 < 1. 

To be more specific, the collision operator / in (JJJ for c ^ 1, acting on exponentially 
bound functions / ~ exp[— /3c fc ], reduces to a linear operator on /, because g ~ Ci and 
c[ ~ ci — p(ci • n)n, yielding the dominant behavior of Ref. [11], and a new important 
subleading correction [26], 

Ioo(c\f) = -f3 a c v [l - /C CT c- b ( CT+1 >/ 2 ]/. (5) 
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Here lC a is non- vanishing, except for 1-D. Its explicit form contains averages {{g s )) and (c s ) 
[26]. The leading term is a direct generalization of the standard methods of Ref. [11] for 
inelastic hard spheres (y = 1) and Maxwell molecules (y = 0). Eq.JSJ allows us to calculate 
for c> 1 sub- leading multiplicative corrections to /(c) of the form, 

In /(c) ~ ~f3c b + [3'c b ' + X lnc + 0(1), (6) 

where b > b' > 0. By inserting |J5J-(JBJ| in 10} , and equating exponents and coefficients of the 
leading and sub-leading terms, we recover the b— values in ©, and determine (3,b',(3',X- We 
only list the values for a = 1 (for further details, see [26]): 

Wwn = Vd(d+l)/( P q((g»+i))) ; (/?6)nf = (d + l)(cf +1 )/(pq({g v+2 ))) (7) 

which includes the special cases considered in Ref. [11]. Regarding the sub-leading quantities 
we find for a — 1, all values of q = ±(1 — a) and d, 

(3' = 0, x = -\v - (d - 1)? 2 /2(1 - q 2 ) (WN) 

/3' = 0, X = -9+(d-l)q 2 /(l~q 2 ) (NF). (8) 

For a > 1 we find sub- leading results similar to (jHJ with xwn = i(l — d) — \v and xnf = 
1 - d- 6. For cr < 1 we obtain b' = ±6(1 - <r), x = 0,0' oc /3/C CT , both for WN and NF 
with corresponding values for b and /3. A striking example of the importance of the new 
sub-leading corrections is shown in Fig^left) for a model with a = —1/2. The numerical 
results are obtained from the Direct Simulation Monte Carlo technique (DSMC), which offers 
an efficient algorithm to solve the nonlinear Boltzmann equation [28]. The dashed-dotted 
curve in Fig. ^left) shows that the simulated c— values are not large enough to agree with the 
asymptotic predictions. However, the solid curve shows a striking agreement with the theory. 
This demonstrates the importance of the sub-leading corrections, which allow to extend the 
validity of the asymptotic theory to much smaller c— value, which are likely to be helpful in 
experimental verification. Such corrections are generic in dissipative gases for d < 2, and more 
pronounced when a < 1. In ID they are absent (see Fig^right)). 

Next we analyze the integral equation (0J at the stability threshold (6 = 0), where power 
law tails, /(c) ~ c~°, are to be expected. The asymptotic approximations, made in are no 
longer valid because of the slow decay of /(C2) in I(c\f) at large velocities. Here we proceed 
as follows. As far as large velocities are concerned (c ^> 1), the thermal range of / may be 
viewed to zeroth approximation as a Dirac delta function <5(c), carrying all the mass of the 
distribution. We therefore write /(c) = 6(c) + h(c), where h represents the high energy tail, 
of negligible weight compared to the thermal bulk. This allows us to linearize the collision 
operator: 

I{c\S + h) = -Ah{c) + 0(h 2 ), (9) 

where we have used that 6(c) makes the nonlinear collision term vanish [Z(c|t5) — 0]. Hence, 
jni defines the linearized collision operator A. Its eigenfunctions are of power law type, and 
particularly suited to describe the algebraic decay of / in cases of marginal stability. For 
the construction of A, which differs from its self-adjoint A^, and of its eigenfunctions and 
eigenvalues we refer to [26]. For the special case of Maxwell molecules the eigenvalues of 
the linearized collision operator have been calculated in Refs. [16,17,27] with the help of the 
Fourier transform method, which can not be applied to non-Maxwell type interactions. The 
essential observation in deriving the spectral properties is that A, applied to any power c a 
generates a term cx c a+u . This property allows us to calculate the eigenvalues X s , and the 
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corresponding left eigenfunctions, A^c s = X s c s+l/ , and right ones, Ac s d v — X s c s d . These 
spectral properties, including 

A S =/3 CT {1-2^1 (-|,^±i;^|l-3 2 )}-^ s+CT , (10) 

are new results. Here 2 F\ is a hypergeometric function, A s is well-defined for s > 0, and 
Ao = is an isolated point of the spectrum. The eigenvalue depends only on the angular 
exponent a, and not on the energy exponent v. Moreover, for s — 2 one has A2 = 2pqf3 G +2 
(see Eq.(@J). Consequently, in the threshold cases (b = 0, which fixes the exponent v at 
threshold to i^wn = — 2 or ^nf = — 1) the asymptotic form of /(c) is given by the right 
eigenfunction that satisfies the integral equation fll|. Substitution yields the transcendental 
equations: 

A s = i } s(s + d-2)X 2 (WN) 

A s = ^sX 2 ({g e+1 ) / (c e+1 ) = ^sX 2 T(9) (NF) (11) 

Their largest root s* determines the exponent of the corresponding high energy tail, /(c) ~ 
l/c s + d + 1 ' for the different driving forces at their stability threshold (6 = 0). 

For the Gaussian thermostat (9 — 1, T(l) = 1) the threshold model (&nf = ^nf = 0) is the 
well-studied Maxwell model with a high energy tail, /(c) ~ l/c s +d , first discovered in Ref. 
[16, 17]. In the 1-D case (where /3 S = 1) the eigenvalue (|10|) simplifies to A s = 1 — q s — p s , and 
(fTTf 1 can be solved analytically for WN, yielding s*(d = 1) = 3 and /(c) ~ l/c s "" +d_2 ~ 1/c 2 , as 
well as for NF (6 = 1), yielding s^ F (d = 1, 6 = 1) = 3 and /(c) - 1/ c s ' '+ d + e ~ 1 ~ l/ c 4 . For all 
remaining (d, ff) values, (|llf> can be solved numerically, where the ratio T(Q) is independently 
measured in the DSMC method. All power law exponents are found to be extremely accurate 
for all driving mechanism, for all values of the restitution coefficient a, and for all dimensions 
(see Fig. 2). So, a new class of power law tails is obtained for WN and NF friction. 

Let us now discuss the power law tail generated by an energy source 'at infinity', as recently 
proposed in Ref. [25] with suggestions for experimental tests. Here the source regularly injects 
a macroscopic amount of energy in the system at ultra high energies. This energy cascades 
down to lower energies, and builds up a stationary 'low' energy tail h(c). The mathematical 
framework developed in the present paper is also well suited to discuss this mechanism. By the 
arguments above I© the collision operator can be linearized, and the tail distribution satisfies, 
Ah = 0. So we determine the right eigenfunction, h ~ c s ~ d ~ v with eigenvalue X s * — 0. As 
X s (s — > + ) = — 1, there does exist a non- vanishing positive solution s*(d, cr, a) > 0. As shown 
in Ref. [26] , the equation for s* is equivalent to the one solved in [25] for the soft sphere model 
(cr = v). In 1-D it has the analytic solution s*(d = 1) = 1 and h(c) ~ c _2_l/ . For d > 1 the 
equation can been solved numerically. Also the more traditional thermostats, as discussed 
in the context of inelastic gases in Ref. [20], can be analyzed with the help of the present 
method [26]. 

In this Letter, we have studied a general class of inelastic soft sphere models, on the basis 
of the nonlinear Boltzmann equation coupled to stochastic or deterministic driving forces. We 
provide a general framework which embeds in a versatile description the majority of driving 
devices and of particle interactions, discussed in the literature: from hard scatterers like 
inelastic hard spheres to soft scatterers like Maxwell molecules; from driving by white noise, 
nonlinear friction, an energy source "at infinity", as well as by traditional thermostats. It 
therefore represents a minimal model for interpreting experiments on driven granular gases, 
discussing trends, etc. When the non-equilibrium steady state is an attractive fixed point 
of the dynamics, the velocity distribution /(c) has a stretched exponential tail oc exp(— c b ); 
sub-leading corrections appear in certain regions of model parameters (v, a, 6), where /(c) is 
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found to be of the form c x exp(— f3v b + (3'v b ). These corrections turn out to be of particular 
relevance since they allow to obtain an agreement with asymptotical theoretical predictions in 
a much larger range of velocities, therefore more easily accessible to experiments. On the other 
hand, algebraic distributions emerge in cases of marginal stability (6 = 0), where power law 
exponents have been calculated. We have argued that such distributions should be extremely 
difficult to observe experimentally. 

Our theoretical starting point is backed up by recent experimental findings, and all predic- 
tions -derived from a new analytical method- have been tested against very extensive Monte 
Carlo simulations [26]. Our work therefore calls for further experimental investigations; in 
particular, electrically forced gases with a charge distribution seem to provide interesting 
test-beds. 
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